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Abstract. We investigate the optical properties of periodic composites containing 
metamaterial inclusions in a normal material matrix. We consider the case where 
these inclusions have sharp corners, and following Hetherington and Thorpe, use 
analytic results to argue that it is then possible to deduce the shape of the corner 
(its included angle) by measurements of the absorptance of such composites when the 
scale size of the inclusions and period cell is much finer than the wavelength. These 
analytic arguments are supported by highly accurate numerical results for the effective 
permittivity function of such composites as a function of the permittivity ratio of 
inclusions to matrix. The results show that this function has a continuous spectral 
component with limits independent of the area fraction of inclusions, and with the 
same limits for both square and staggered square arrays. 
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1. Introduction 

This paper links themes evoked in two classic papers, one in mathematics pj| and the 
other in physics [2]. The first of these poses the question as to whether the spectral 
content of the radiation from a body can reveal its shape. The second shows that the 
use of negatively-refracting metamaterials in a plane slab can lead to a super-resolving 
perfect lens, also known as a superlens. We will consider a two-dimensional composite 
material, composed of polygonal inclusions made of a metamaterial (by which we mean 
an artificial material with a dielectric contant which has a negative real part and a 
very small imaginary part) and placed in a positive dielectric matrix material. We will 
show that, in the spirit of Pendry, the metamaterial makes possible resolution of an 
important structural feature of the inclusions, irrespective of how much smaller than 
the wavelength they are. We will also show that, in the spirit of Kac, this feature relates 
to the shape of the inclusion, being in fact the corner angle of the polygon, and that 
it is deduced from spectral measurements on the composite. The fact that a spectral 
feature could be determined by corner shape, independent of (say) the area fraction 
of inclusions, was first suggested by Hetherington and Thorpe on the basis of an 
elegant argument and numerical evidence for dilute composites. 

We will base our demonstration firstly on analytic results relating to the spectrum 
of the effective dielectric permittivity function ees of the composite material, and 
secondly on remarkably accurate numerical results for this spectrum obtained using 
a new technique. Note here that we are using the word "spectrum" in two related, but 
slightly different senses. In the previous paragraph, its usage meant that the absorption 
of electromagnetic waves by the composite was being determined as a function of 
wavelength, ranging over an appropriately-wide band. In the first sentence of this 
paragraph we referred to the distribution of singularities of the function eefr, giving 
the effective permittivity of a composite having a specified geometry as a function of 
the ratio a = ei/e2 of the permittivities of inclusions and matrix. The relation between 
these usages is that, as wavelength varies so does the ratio a, so that measurements of 
(say) optical absorption by a composite over a suitable wavelength interval can reveal 
details of the singularities of the function egfr. 

The numerical results for the singularity spectrum of eefr reveal that it has a 
continuous part which runs between upper and lower limits of a which do not vary at 
all with the area fraction of the inclusions. It is complemented by a discrete spectrum 
of poles which does evolve with area fraction. This evolution is in fact necessary, since 
the continuous spectrum for touching square inclusions in a checkerboard arrangement 
occupies the entire negative real axis of a, but for non-touching square inclusions is 
confined to the interval —3 < a < —1/3. The animations we give show how this 
transformation is achieved: the discrete spectrum becomes more and more dense as 
the touching configuration is approached, to supply the required spectral extension, as 
anticipated by one of the authors 

The results we give here are interesting in the insights they give into the connection 
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between metamaterials and super-resolution. They are also important in furthering 
our understanding of the connection between inclusion shape, geometrical arrangement 
and spectral properties of the effective permittivity function. This connection helps in 
the design of structures having enhanced absorption over a wide wavelength range for 
applications in photothermal or photovoltaic captors [SI El El E], or offering strongly 
enhanced local fields for applications like sensing or nonlinear optical elements. 

We give a brief overview in Section |2] of some of the important properties of the 
function eefr(o"), including analytic results relating to the continuous part of the spectrum 
and some numerical investigations of both the discrete and continuous parts of the 
spectrum. In Section [3} we describe the method which enables accurate calculation of 
the spectrum of ees, and give numerical results illustrating its convergence, even on the 
negative real axis of a. In section |4], we discuss the animations which are given in the 
Supplementary Material to this paper, and the physical consequences of the behaviour 
they show. We give a discussion and concluding remarks in Section [5j 

2. Properties of the effective permittivity function for composites 

We will now give a concise review of what is known about the properties of the effective 
dielectric permittivity function ees, for composites made of two materials with dielectric 
permittivities ei and 62, with the former corresponding to a disconnected inclusion 
phase and the latter to a continuous matrix phase. We assume the geometry has cubic 
symmetry (in three dimensions) or square symmetry (in two dimensions) so that is 
scalar valued, i.e. the effective dielectric tensor equals eeg/. This review builds on that 
given in Perrins and McPhedran [9]. 

The calculation of eefr for a given geometry is homogeneous of degree 1 in the 
variables ei and €2, so we can rescale to make eefr a function of a single complex variable, 
the permittivity ratio cr = exje^. 



From now on we let eefr be an abbreviated notation for the effective dielectric permittivity 
function eefr(ei,e2) and we let eefr(cr) denote its scaled counterpart, the effective relative 
dielectric permittivity function, as defined via ([T]). 

To determine eefr(cr) for a given geometry, we need to solve an electrostatic transport 
problem repeatedly for various cr. More precisely, we need to solve Laplace's equation 
for the potential on a periodic domain with a periodic electric field — VV^(x, y) having 
a prescribed average value, and boundary conditions of continuity of V and its normal 
flux edVdn at interfaces between materials. The theory of the function eefr(a") becomes 
particularly elegant when we deal with two-dimensional problems, in which V = V{x, y) 
becomes a function in the plane. We may then apply the apparatus of complex-variable 
theory to the calculation of V , and thus to ecfT(cr). For the case of a doubly-periodic 
array of inclusions C with unit cell W, and square symmetry, the effective permittivity 




(1) 
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may be defined as [10 



\J^\/Vix,y)dxdy\^' 

where the integral in the numerator includes contributions Si from the inclusion region 
and S2 from the matrix region. Except for occasional comments on effective permittivity 
in three-dimensions, we will concentrate on two dimensions, which corresponds to arrays 
of cylinders of arbitrary cross-section, with the average field aligned in the (x, y) plane. 
The area fractions of the two phases will be denoted pi and p2- 

Since the geometry has square symmetry, Keller's Theorem [11] gives 

eefr(cr)ecfr(l/o-) = 1 • (3) 

This equation then pairs zeros of eeg(cr) at values ctq with poles at values cXp = I/cq- 
Of course, from ([s]), zeros of ecfr(o") require that the contributions £1 and £2 add up 
to zero. Since, with €2 = 1, £2 is real and positive, this means £1 must be real and 
negative, and so €2 = cr must be real and negative at any zero of eee{<7), and thus 
at any pole as well. Bergman [10] proved this property for both two-dimensional and 
three-dimensional composites. Bergman also recognized that even though the transport 
problem when the average electric field — VV^(x, y) is prescribed does not have a solution 
at a pole, it should have a solution at a pole when instead the average displacement 
field — eVV^(a;,?/) is prescribed. Milton [12] proved that for the value a = —1 the 
electrostatic problem of an array of circular cylinders, with either a prescribed value of 
the average electric field — or a prescribed average value of the displacement 

field —eWi^Xjy), does not have a solution, compare with the discussion of (A.l) below. 
This suggests that a = — 1 is either an essential singularity or lies on a branch-cut of 
the function eefr(cr) for arrays of circular cylinders. 

The fact that branch-cuts cannot be in the upper or lower half-planes, but must lie 
exactly on the negative real axis of a was proved by one of the authors [13|, using the 
relationship between composite materials and resistor networks. A rigorous justification 
of the spectral representation for eefr(o-) was given by Golden and Papanicolaou [H]. In 
general the function ecs{<y) has the representation 

eefr(a) = ao + aicr+ / — (4) 

where Oi and the spectral measure dfi{T) are non-negative. The support of dfi{T) is the 
spectrum. The spectral measure can be recovered from the values that the imaginary 
part of eefr(o") takes near the negative real cr-axis since the integral of any smooth test 
function g^r) with respect to the measure (i/i(r) is given by 
.0 ^ .0 

/ g{T)dn{r) = hm -/ g{T)Qecs{r + 15) dr. (5) 
J-00 5-^0^ J -00 

The discrete spectrum of eeff(cr) is readily exhibited numerically. This has been done 
for arrays of spheres by Bergman [TJ] and for arrays of circular cylinders by McPhedran 
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and McKenzie [16], with both studies showing that the poles and zeros of eoff(cr) converge 
to an essential singularity at a = —1. 

We focus now on what can be said about the continuous spectrum of ees{cr). One 
simple geometry for which a result is immediately apparent is the checkerboard, for 
which Dykhne [T7] obtained from Keller's theorem ([s]) the exact result 

eefr(cr) = A/a. (6) 

This then exhibits a branch cut along the entire negative real axis of a. 

An exact result can also be obtained for the polarizability a of a pair of touching 
cylinders [5]. Using an inversion of coordinates about the contact point, the touching 
cylinders may be transformed into a slab of matrix material with permittivity €2 = 1 
surrounded by two half-planes filled with material with permittivity ei = a. Introducing 
the parameter 

A = ^, (7) 

o- + 1 

it is easy to show using the method of images that the polarizability for a pair of touching 
cylinders of radius a for the case of the applied field parallel to the line connecting 
cylinder centres is 

°° A' 

a = 4^2^-. (8) 

1=1 



We see that the series in ([8j) converges provided |A| < 1, i.e. for real a, a > 0. However, 
we can obtain a meaningful result even when this is not the case by the technique of 
analytic continuation, since the series in ([s]) is a known transcendental function, called 
the dilogarithm, and denoted Li2. Thus, we can replace ([s]) by 

a = 47ra2Li2(A) . (9) 

The properties of the dilogarithm function are that it has a branch cut running from 
A = 1 to A = oo, across which the discontinuity in the imaginary part of Li2(A) is 
27r log[3?{A}]. The branch cut in the plane of relative permittivity runs from a = —oo 
to 0" = — 1. If the direction of the applied field is perpendicular to the line of centres, 
the branch cut runs from o" = — ltoo" = 0. 

We next consider arrays of square inclusions, for which we have already mentioned 
the Dykhne result ([6]). A generalization of this for an array in which the square unit 
cell was divided into four equal squares with dielectric permittivities ei, €2, 63, and €4, 
was conjectured by Mortola and Steffe [18]: 

Ceff = 



1 /2 / \ 1 /2 

(e2 + e3)(e4 + ei)] / eie2e3 + £16264 + eie3e4 + e2e3e4\ 



/ eie2e3 + ^1^2^ + ei^sQ + e2e3e4 \ ^^^-^ 
(ei + e2)(e3 + €4)] V ei + €2 + €3 + 64 / 

This conjecture was proved independently by Milton [ID] and Craster and Obnosov f2U\ . 
We will consider a particular sub-case of this result, due to Obnosov [2T]: an array 



of square cylinders with area fractions pi = 0.25 and p2 = 0.75, for which (10) gives 



, , ; 1 + 3cr ^ ^ 

eef^(a) = ^/^^. (11) 
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This formula yields a spectrum consisting solely of a branch-cut running from a = —3 
to a = —1/3. We compare the result given by this formula for the real part of eefr(cr) 
with the result of a numerical mode matching procedure in Fig. [1} This comparison 
reveals the difficulty of evincing details of the spectrum using numerical methods: the 
mode matching technique approximates the branch cut by a discrete set of poles, which 
becomes more dense as the number of modes increases. However, it is difficult to 
distinguish between branch cuts and sets of poles concentrating around an essential 
singularity by such methods. Furthermore, the mode matching method failed to give 
clear indications of the spectrum for area fractions of cylinders distinct from pi = 0.25. 

The behaviour of fields near corners of inclusions in a matrix of differing 
dielectric permittivity and magnetic permeability was treated in electromagnetism by J. 
Meixner [22]. As well as obtaining the exponents which characterize the singularity of 
the field components near the edge, Meixner pointed out that exactly the same formulae 
could be applied in electrostatics and magnetostatics. Hetherington and Thorpe [3] 
analysed the behaviour of fields near corners in electrostatics and magnetostatics, 
apparently without knowledge of Meixner's paper. They however pointed out the link 
between field behaviour near corners and the nature of the singularity spectrum in 
composites containing inclusions with corners. Let us suppose that the electrostatic 
potential varies with distance r from a corner with an included angle 2ip as r^, and that 
the permittivity ratio inside the inclusion to that outside is a. Then [221 [3] /3 is found 
by solving a transcendental equation: 
tan[(3 {tt - iP)] _ 

tan(/3V^) " ^ ' 
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or the same equation with a replaced by l/a. For a real, the solution /3 of ( 12 ) is either 
pure real or pure imaginary. 

To clarify the behaviour of fields, consider the case oi i{j = 7r/4, corresponding to a 
90° corner. Then (12) gives 

From this, we see that the right-hand side exceeds one in magnitude for a lying between 
—3 and —1/3, and that /3 is then pure imaginary, corresponding to a solution for the 
potential which oscillates with r, and the oscillations become ever more rapid as r tends 
to zero. The electric field is given by the spatial derivative of the potential, and it 
diverges like 1/r multiplied by the oscillating term as r — > 0. The branch cut region 
here is that where (3 is imaginary, i.e., between a = —3 and cr = —1/3. Hetherington 
and Thorpe [3] postulate that for an m sided regular polygon, there will be a branch 
cut running between a = (2 + m)/(2 — m) and a = (2 — m)/(2 + m). They made this 
assertion after recognizing that in this interval of a the surface charge (rather that just 
the surface charge density) near the corner is infinite, which is unphysical. 

Another argument for the position of this branch cut was put forward by one of 
the authors jlj. In order that ees have a significant imaginary part when e has a very 
small imaginary part we see from ^ that the electric field — must be close to losing 
its square integrability. From the asymptotic form of V near the corner one finds that 
this happens when the imaginary part of a is very small and real part of a is between 
(2 + m)/(2 — m) and (2 — m)/(2 + m). To elucidate this further for the case of a 90° 
corner we solve ( 13 ) with P = /3' + and a = a' + ia" where /3', a' and a" are real, 
—3 < cr' < —1/3, and a" and /3' are both very small and positive. This gives 

4a" 

(5' ^ ^^^^^=. (14) 

((T' + l)v/(a' + 3)(3a' + l) 

Using polar coordinates (r, 6) near the corner, the potential V scales as r'^ as r — )■ 0, and 
so |Vl^p will be close to r^^'~^|(7(6', cr')p for some function g{6, a'). It follows that with 
ei = 0" and €2 = 1 the imaginary part of Si (which is a measure of the power dissipation 
in the composite) has a contribution near the corner from inside the radius r = rg of 

de I a"\VV\^r dr = j \g{e,a')\^ dO 

213' 

r, 



^|(a' + l)v/(a' + 3)(3a' + l)| J \g{e,a')f dO 

(15) 



where the integral over 6 is only over those angles in the inclusion. Thus, provided 
g{6, cr') is non-zero, the contribution of the corner to the imaginary part of Si remains 
non-zero even in the limit a" — t- (and goes to zero when a' < —3 or cr' > —1/3 since 
then a"//3' 0). 

A corner is not the only geometric feature which can act gnificant energy 

absorber when the imaginary part of the dielectric constant goes to zero. The center of 
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Figure 2. Left: a cutout from C of a square array of square cylinders witli area 
fraction pi ~ 0.5 and a unit cell U with a 16-panel coarse mesh on Fg. Right: the 
same thing for a staggered array of square cylinders, but with pi — 0.4 and a 32-panel 
coarse mesh on Tq. 

a sphere with a dielectric constant Ai in the radial direction and dielectric constant A2 
in the tangential direction acts as an absorber when A2/A1 approaches real values less 
than —1/8: see figure 4 in the paper of Qui and Lukyanchuk [24j (which shows that this 
energy absorbing feature extends beyond the quasistatic limit) and see also the related 
discussion on page 239 of |1] (where there is an error as the limit 6 — )■ —1/2 should have 
been taken rather than the limit 5 — )■ 0). 

3. Numerical method 

We now describe a numerical method stable and accurate enough to verify the conjecture 
of Hetherington and Thorpe for the case m = 4, and to show the spectral evolution as 
a function of area fraction for composites with square inclusions. 

Laplace's equation is to be solved on a doubly periodic domain C. The boundary 
conditions on the positively oriented interface F between the inclusion phase and the 
matrix phase are given in Section [2j An average electric field Eq of unit strength 
is applied. The permittivity of the matrix phase is set to €2 = 1 so that the effective 
permittivity is equal to the effective relative permittivity. From the repeated solution to 
this problem for various ei = cr we obtain eeff(o"). Two types of domains are investigated: 
the "square array of square cylinders" and the "staggered array of square cylinders" (a 
square array of cylinders with diamond-shaped cross-sections, as studied for example in 
[251 EH]), see Fig. g 

Boundary value problems on domains involving sharp corners may require extreme 
resolution close to corner vertices, even when the demands for overall accuracy are 
moderate. One has to be selective with the choice of numerical method. We chose an 
integral equation based scheme. Such schemes have the advantage that they can retain 
stability also in very difficult situations. 

Our particular choice of integral equation is standard - a single-layer equation [27] . 
For its solution we use a novel numerical method called recursive compressed inverse 
preconditioning. Conceptually this is a local multilevel technique which makes a change 
of basis and expresses the non-smooth solution to the single-layer equation in terms 
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Figure 3. Left: the effective relative permittivity of a square array of square cylinders 
at pi = 0.25. The curves are supported by 826 adaptively spaced data points (not all 
values shown due to the setting of the axes). Right: the relative error with (11) as 
reference value. 



of a piecewise smooth transformed layer density which can be cheaply resolved by 
polynomials. Discretization leads to a block diagonal transformation matrix R (an 
inverse preconditioner) where the columns of a particular block can be interpreted as 
special basis functions for the original density in the vicinity of a corner vertex multiplied 
with suitable quadrature weights. The blocks of R are constructed in a fast recursion, 
i = l,...,n, where step i inverts and compresses contributions to R involving the 
outermost quadrature panels on level z of a locally n-ply refined mesh. We emphasize 
that the method is strictly numerical and fully automatic. There is no separation of 
variables or eigenvalue analysis involved. 

The recursive compressed inverse preconditioning method was originally described 



in Ref. [28] and further developed in Refs. [291 EOl E]. Appendix A below highlights 
some of the method's features, relevant to the domain C. A fuller description will be 
included in a forthcoming paper |32j . 

3.1. Achievable accuracy 

We first compute eeff(cr) for the square array of squares at pi = 0.25 in the limit of a 
approaching the negative real axis from the upper half-plane H, as in the example of 



Fig. 1. The exact result (11) is used as a benchmark. Fig. |3j shows that the relative 
error is close to machine epsilon (the upper bound due to rounding in floating point 
arithmetic) except for in a neighbourhood of three points where it is higher: the ends of 
the branch cuts at cr = —3 and a = —1/3, and at the singularity of the integral equation 
at cr = —1. This demonstrates that the problem of computing eeff(cr) for arrays of square 
cylinders is well conditioned in general and that our scheme is stable. 

The staggered array of square cylinders at pi = 0.49999999995 is a more challenging 
geometry than the square array of square cylinders at pi = 0.25: 
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Figure 4. Left: the effective relative permittivity of a staggered array of square 
cylinders at pi — 0.49999999995. The curves are supported by 2006 adaptively spaced 
data points. Right: the absolute difference between the left- and the right hand side 
of (fiel). 



• there are more length scales involved, 

• eefr(c") varies more rapidly and has more poles and zeros, 

• there is no exact result to compare with. 

The first problem is the least difficult. The multilevel property of our numerical method 
should enable the resolution of almost arbitrarily small separation distances between 
corner vertices. The second problem is more severe. For eeff(o') close to zero, one can 



expect cancellation in (A.2) and the relative accuracy should suffer. Furthermore, it 
is harder to resolve wildly varying functions in floating point arithmetic than slowly 
varying ones. The third problem is solved by using the extent to which ([s]) is satisfied 
as an indicator of the relative error. For this, since a and 1/a lie on different sides of 
the real axis and our numerical method takes limits from H, we use ([s]) in the equivalent 
form 

eefr(a)e:ff(l/a*) = l, (16) 

where the symbol to denotes complex conjugation. Fig. |4] suggests that despite the 
difficulties, typically, only a few digits are lost compared to the square array at pi = 0.25. 
The numerics seem to give a relative precision of at least 10~^ even for the most extreme 
values of eefr(cr). 

For the interpretation of various limits it might also be of interest to study ecfr(o") 
for a some finite distance into H. Fig. [5] shows again the staggered array of square 
cylinders at pi = 0.49999999995, but unlike in Fig. |4]we have here interrupted the limit 
process at a a relative distance of 10"^ away from the real axis. 
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Figure 5. Same as Fig. |4j but a is multiplied with a complex constant 1 — i • 10^^. 
The curves are supported by 3435 data points. The red curve is drawn on top of the 
blue curve. 



4. Animations of Spectral Evolution 

This section discusses tlie evolution of the variation of ees{<^) with permittivity 
ratio a for the square array of square cylinders and the staggered array of square 
cylinders as the area fraction of squares, pi, varies. The discussion relates 
to two animations, called Animation 1 and Animation 2, which can be viewed 



at |http : //www .maths . 1th . se/na/ staff /helsing/animat ions . html To facilitate 



viewing, each animation is available in four versions, denoted A, B, C, and D. The versions 
have the same content, but differ in image format and pixel resolution. 

We first consider the evolution of the variation of eeff(cr) as pi in the square array 
ranges from zero to unity. This evolution is shown in Animation 1, from which a typical 
frame is given in Fig. |6] The most important feature of Animation 1 is quite clear: for 
a real, and for all values of pi, and except at the poles, non-zero values of 3{eefr (cr)} are 
confined to the interval —3 < a < —1/3, in accord with the suggestion of Hetherington 
and Thorpe [3J. Of course, the value of this imaginary part is always positive, if we 
restrict ourselves to composites without gain (for which the imaginary part would be 
always negative). 

Below the area fraction of 0.25, the real part of eefr(cr) is positive, and it develops 
its first pole at this area fraction. It is interesting to compare frames from Animation 1 
and the left image of Fig. [3] with Fig. [TJ the results of the mode-matching method clearly 
correspond to those of the new method, but are capable of a resolution limited by the 
number of terms employed in field expansions. 

For pi > 0.25, the real part is negative between a = —3 and a = —1/3, while 
the pole migrates to more negative values of cr. For area fractions near pi = 0.75 (see 
Fig. [6]), "features" which we call quasipoles develop from near a = — 1, and one moves 
towards cr = —3, while the other moves towards a = —1/3. When they reach these 




Figure 6. Real (blue) and imaginary (red) parts of the effective relative dielectric 
permittivity for a square array of square prisms with area fraction pi = 0.735 as a 
function of permittivity ratio a. 



values, and then are not muted by the absorbing nature of the corners, they transform 
into actual poles, which move towards a = — oo and a = respectively. At higher 
values of area fraction, more quasipoles evolve from a = —1 and give rise to additional 
actual poles when they emerge from the branch-cut region. 9{eeff(cr)} becomes small 
as pi — )■ 1, while 3fJ{ecff(o")} tends towards a, apart from the increasingly numerous but 
increasingly narrow pole regions. 

For the staggered array of square prisms. Animation 2 illustrates the behaviour of 
eefr(cr) as a function of a for area fractions ranging from zero to 0.5 (with the behaviour 
for pi in the range 0.5 to 1.0 following from that in the lower range using Keller's 
Theorem ([3]). As has been commented in Section [2| the interesting question is how 
the branch-cut location (from cr = — oo to cr = 0) for pi = 1/2 of equation ^ can be 
reconciled with that (from o" = — 3tocr = — 1/3) predicted by Hetherington and Thorpe 
[3j for pi arbitrarily near 1/2. A mechanism for this reconciliation was provided by one 
of the present authors [1]: discrete sets of poles in — oo < a < — 3 and 1/3 < a < were 
predicted to become denser and denser as pi approached 1/2, thus extending the branch 
cut in the limit to that required by ([6]). The accuracy of this prediction is evident in 
Animation 2: poles develop from the quasipoles generated at a = —1, and move left 
and right into the embryonic branch-cut regions — oo < <j < —3 and 1/3 < o" < 0. 
The left frame in Fig. 4 shows a stage in this evolution where pi is very close to 1/2. 
Animation 2 makes the "nursery role" of the region around a = — 1 in the development 
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of the spectrum much more evident (due to larger amphtudes of the quasipoles) than 
does Animation 1. 

The sensitivity of the spectral details for the staggered array near the checkerboard 
configuration are very evident in Animation 2. As we have commented in Section |2| the 
asymptotics of fields near corners are the same in electromagnetism as in electrostatics 
[22]. Thus, attempts such as that in |23] to model the transition from electromagnetically 
reflecting structures to electromagnetically transmitting structures as pi moves through 
1/2 would require an adaptive and recursive method hke that described in Section [s] to 
be able to achieve sufficient accuracy. 

5. Discussion and Conclusions 

In this paper we have brought together rigorous mathematical results with numerical 
investigations of unprecedented accuracy. The latter have revealed the generality of the 
former, and have substantiated a conjecture of Hetherington and Thorpe [3^ in a striking 
and conclusive way. 

We conclude by commenting further on how the arguments and results we have 
presented can be implemented in a practical demonstration of morphological super- 
resolution, uniting the ideas of Kac [1] and Pendry [2]. Such a demonstration would 
require the fabrication of a set of parallel cylinders with a square cross section (or 
polygonal cross section). The cylinders do not have to be arranged in a geometrically- 
perfect array, and they do not have to be densely packed. They have to be made 
of a material which is essentially non-absorbing and with a negative permittivity or 
permeability. 

These requirements suggest the set of cylinders be made of metal, of size 10 /im 
or larger, and be probed with wavelengths far greater than the cylinder size in the far 
infrared or longer. Such cylinders are large by today's lithographic standards, and so it 
should be possible to accurately form their corners to achieve sub-wavelength accuracy. 
Going into the far infrared region diminishes metallic loss from its value in the visible and 
near-infrared [55]. It is crucial that the metallic loss be very low, since the experimental 
signature we suggest be probed is enhanced absorption by a set of such cylinders over a 
wavelength interval in which the metal's permittivity ranges from say -1/3 to -3 (scaled 
relative to the permittivity of a host dielectric in which the cylinders are embedded). 
Note that the enhanced absorption of incident radiation detected will increase as more 
lines of cylinders are added to the set. 

The experimental result which would indicate morphological super-resolution is an 
enhanced absorption for wavelengths far in excess of the cylinder size, switching on and 
off at geometrically-determined limits described above, independent of the arrangement 
and area fraction of the cylinders. We stress however that such a demonstration would 
indicate the physical relevance of the ideas we have described for a particular system. 
The mathematical results we have described are of course rigorous, and the numerical 
examples of them we have given are highly accurate, so our demonstration of super- 
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resolution for wavelengths arbitrarily larger than the size of the particles probed does 
not rely for its validity on experimental support. They may be applicable to governing 
equations other than the Helmholtz equation, for which the ideas of metamaterials and 
their applications are currently being explored [36] . 



Acknowledgments 

Ross McPhedran acknowledges support from the Australian Research Council's 
Discovery Projects and Centre of Excellence Schemes. Graeme Milton acknowledges 
support from the National Science Foundation through grant DMS-0707978. 



Appendix A. Key features of the numerical method 

To keep the notation short we make no distinction between points or vectors in a real 
plane and points in a complex plane C. All points will be denoted z or r. 



The integral equation 



The potential function V{z) in C is represented as a sum of a driving term and a single- 
layer potential with density p{z) [33]. Enforcement of the boundary conditions on F 
leads to the Fredholm second kind integral equation 

\ f ( n n* dr 1 

p(z) + -lp(r)'^\^^\=2X^{E;n,}, zeTo- (A.l) 



Here is the outward unit normal of F at z, Fq denotes the restriction of F to U, 
and A is as in ([t]). Note that, as a — )■ — 1 we have A ±oo and (A.l) is no longer 



a second kind equation, but a first kind equation whose (unique) solvability is by no 
means guaranteed. Therefore one can say that a = —1 is a singularity of (A.l). 



Once (A.l) is solved for p(z) and under the assumption that the inclusions do not 



overlap the unit cell boundary, the effective relative permittivity in the direction of the 
applied electric field can be computed from 



eefT(Cr) 



p{z)^{E;z} d\z\ . 



(A.2) 



ro 



Depending on how the unit cell is chosen, the squares in the staggered array may 
overlap the unit cell boundary. With the choice in Fig. [2| they certainly do. But since 
p{z) is a periodic function and identical on all squares one can circumvent this problem 



by modifying (A.2) so that it integrates p{z) twice on the square at the center of the 



unit cell and ignores p{z) on the other squares. 



Discretization 



We discretize (A.l) and (A.2) using a Nystrom method based on composite 16-point 
polynomial interpolatory quadrature and a parametrization z{t) of F. The parameter t 
is real. See Ref. |35l fo^' ^ review of Nystrom methods including error analysis. 
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An initial coarse mesh that resolves the kernel of the integral operator in (A.l) 



away from the corner vertices is constructed on F, see Fig. [2j The coarse mesh is refined 
by subdividing those panels that neighbour corner vertices. The subdivision is done 
n times in a direction towards the vertices. On quadrature panels which neighbour 
corner vertices we choose quadrature nodes according to the zeros of certain Jacobi 
polynomials. On remaining panels we choose quadrature nodes according to the zeros 



of Legendre polynomials. Upon discretization on the refined mesh (A.l) assumes the 
form 

(Ifine + Kfine) Pfine = gfine , (^-3) 

where Igne and Kgne are square matrices and pgne and ggne are column vectors. The 
vector gfine corresponds to the discretization of the piecewise smooth right hand side. 



Now the kernel K{t, z) of the integral operator in (A.l) is split into two functions 

K{t,z) = K\r,z) + K°{T,z), (A.4) 
where K*(t, z) takes care of corner interaction and K°(t, z) can be viewed as the kernel 



of a compact integral operator. The kernel split (A.4) corresponds to an operator split 
and the change of variables 

p{z) = {I + KT' K^) (A.5) 



makes (A. 3) read 



(ifine + Kfi^g (Ifine + ^hnc) ) P& 



gfir 



(A.6) 



This right-preconditioned equation corresponds to the discretization of a Fredholm 
second kind equation with compact operators. The solution pfine is the discretization of 
a piecewise smooth function. 



Compression 



The matrix Kg^^, the density Pfine, and the right hand side gfine in (A.6) can be evaluated 



on the coarse mesh without the loss of precision. Only (Ifine + Kg^^) ^ needs the refined 



mesh for its accurate evaluation. This enables a compression of (A.6). We introduce 
the compressed weighted inverse 



R 



Ifine + K 



,-1 



fine/ 



P. 



(A.7) 



Here P is an unweighted prolongation operator that performs panelwise 15th-degree 
polynomial interpolation in the parameter t (which as we recall parameterizes F through 
z{t)) from points on the coarse mesh to points on the fine mesh when acting on column 
vectors from the left. Pw is a weighted prolongation operator. See Section 5 of Ref. [30] . 



Substitution of (A.7) into (A.6) and the use of some relations between prolongation 



operators make (A.6) assume the form 



(Icoarse + -Kcoarse-f^) Pc 



(A. 



This equation, defined solely on the coarse mesh, will be used in our computations. 
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Fast recursion for R 



The construction of R from its definition (A. 7) may be a costly and unstable operation 



when the refined mesh has many panels. The number of subdivisions n needed to reach 
a given accuracy may grow without bounds due to the singularities in p{z) that arise as 
cr approaches certain values. 

Fortunately, the construction of each block of R, associated with a corner of the 
square array or with a corner-meet of the staggered array, can be greatly sped up and 
also stabilized via a recursion. This recursion uses matrices K on local meshes centered 
around corners or corner-meets. It would be going too far to describe all the fine details 
of this procedure, but Sections 3.2 and 3.3 of Ref. [31] give a fairly good idea of how the 
recursion is set up in the present context. A key step is the (partial) conversion of the 
recursion into a non-linear matrix equation. This equation is solved using a variant of 
Newton's method relying on numerical homotopy to approach purely negative a from 
the upper half-plane H. The ratio a, which enters into K, is initially multiplied with a 
constant q = 1 — O.Oli. The imaginary part of q is reduced with a factor of ten after each 
of the first 14 Newton iterations. Then q is set to unity and the iterations are continued 
until either a sharp convergence criterion is met or a total of 30 iterations is reached. A 
full description is given in Ref. [32j . 
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